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We use determinantal quantum Monte Carlo simulations and numerical linked-cluster expansions 
to study thermodynamic properties and short-range spin correlations of fermions in the honeycomb 
lattice. We find that, at half filling and finite temperatures, nearest-neighbor spin correlations can 
be stronger in this lattice than in the square lattice, even in regimes where the ground state in 
the former is a semimetal or a spin liquid. The honeycomb lattice also exhibits a more pronounced 
anomalous region in the double occupancy that leads to stronger adiabatic cooling than in the square 
lattice. We discuss the implications of these findings for optical lattice experiments. 
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In recent years, the isolation of graphene flakes |l| 
has generated a revolution in solid state physics |2|. 
Graphene is an atom thick structure with carbon atoms 
arranged in a honeycomb lattice geometry, which features 
low energy excitations that are massless Dirac fermions. 
Given its reduced coordination number, graphene has 
also opened a new venue to create exotic quantum phases. 
Based on quantum Monte Carlo (QMC) simulations of 
the half filled one-band Hubbard model, Meng et al. Q 
have argued that a spin-liquid ground state may be real- 
ized in this lattice geometry at intermediate interaction 
strengths. Earlier works had found the ground state to be 
a semimetal in the weakly-interacting regime and a Mott 
insulator with long-range antifcrromagnetic (AF) correla- 
tions in the strongly-interacting regime 0, HI ■ The find- 
ing of an intermediate spin-liquid phase, recently chal- 
lenged by another QMC study that considered larger lat- 
tice sizes [p| , motivated much theoretical work in related 
models ^UM- 

Experiments on fully tunable artificial graphenc-like 
lattices now offer a pathway to study the physics above, 
and more, in a controlled way (l6Hl9j . Motivated by 
those experimental advances, especially by the availabil- 
ity of an ultracold lattice fermion setup [18| , where on-site 
interactions, hopping amplitudes, doping, and tempera- 
ture can be fully controlled using Feshbach resonances, 
changing the lattice depth and the number of fermions in 
the gas, and varying the cooling time [2(|, respectively, 
we study thermodynamic properties and short-range cor- 
relations of two-component correlated fermions in the 
honeycomb lattice. 

We show that such a system exhibits several unex- 
pected properties when compared with its square lattice 
counterpart. For example, at a finite temperature, it 
may be less compressible in the weakly interacting regime 
where its ground state is a semimetal, while the latter 
becomes less compressible in the presence of strong in- 
teractions when both lattices have an insulating ground 



state. We also identify temperature regimes in which, 
surprisingly, (i) nearest-neighbor (NN) spin correlations 
are stronger and (ii) a more significant anomalous region 
can be seen in the derivative of the double occupancy 
with respect to temperature, in the honeycomb lattice 
than in the square lattice. 
We consider the one-band Hubbard Hamiltonian 
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t ^2 ^\<jC ja + H.c.) + U^2 "*t"*4. . 
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where standard notation has been used [2l|. At half fill- 
ing, in the square lattice, the ground state of this model 
is an AF Mott insulator for any U > [2l[ , while, in the 
honeycomb lattice, it has been recently argued to be a 
semimetal for < U/t < 3.5, an AF Mott insulator for 
U/t > 4.3, and a gapped spin liquid in between 

In this Letter, to study the properties of the Hamilto- 
nian [Eq.([T])] in the honeycomb and square lattices, we 
utilize two unbiased computational approaches, the de- 
terminantal quantum Monte Carlo (DQMC) technique 
-l24j and numerical linked-cluster expansions (NLCEs) 
\2l\ . DQMC simulations are performed in finite-size 
systems (with 100 and 96 sites for the square and hon- 
eycomb lattices, respectively) using a small discretized 
imaginary time (At x t — 0.05). NLCE calculations, on 
the other hand, provide exact results in the thermody- 
namic limit but converge down to a temperature that 
is determined by the divergence of correlations and the 
largest cluster sizes that we can consider. Here, we in- 
clude clusters up to the ninth order in the site expan- 
sion and use Wynn and Euler resummation algorithms to 
extend the region of convergence to lower temperatures 
[M 0, HI- DQMC calculations and NLCEs are com- 
plementary as the former provides more accurate results 
down to lower T for U < w, where w is the nonintcr- 
acting bandwidth (w = 6t for the honeycomb lattice and 
w = 8t for the square lattice) while the latter is better 
suited for U > w [2§|. In the region where DQMC sta- 
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tistical errors are small and NLCEs converge, we obtain 
an excellent agreement between both approaches. 

In optical lattice experiments, single site addressability 
30l 31 1 makes possible an accurate determination of the 
equation of state [density (n) vs chemical potential (/x)] 
of lattice Hamiltonians of interest. This equation of state 
determines the shape of the experimental density profiles 
and, when obtained at low enough temperatures, allows 
one to identify the presence of a single particle gap in 
the spectrum. In the inset of Fig. HJa), we show the 
equation of state in the square and honeycomb lattices 
for U/w = 3/2, which is beyond the critical value for the 
formation of the Mott insulator in the latter, and for two 
values oiT/w that are very close in both lattices. With 
decreasing temperature, nvs/j reveals the single-particle 
gap in the Mott phase by exhibiting a region in which 
n barely changes when changing (i. As expected from 
their phase diagrams, that gap is greater in the square 
lattice than in the honeycomb lattice. This results in the 
former system being less compressible than the latter at 
half filling and finite T for large values of U. 

By decreasing T for small U, the compressibility (k = 
dn/d/j.) also reveals the vanishing of the density of states 
in the semimctallic phase. This is shown in the main 
panel of Fig.[TJa), where, for weak interactions, the com- 
pressibility in the honeycomb lattice is seen to decrease 
with decreasing temperature (k — >• as T — > 0). This 
behavior is to be contrasted with the one in the square 
lattice, where k increases as U — >• and T — > 0, signaling 
the metal insulator transition 2l|. Note that, for finite 
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T, the behavior above leads to a region in U where k 
is smaller in the honeycomb lattice than in the square 
lattice, despite the fact that in such a region the ground 
state in the former may be a semimetal while in the lat- 
ter is an insulator. This can be understood given the dif- 
ference between dispersion relations in the two systems 
which, at low T, can lead to less states being available in 
the honeycomb lattice than in the square lattice. 

Another quantity of much interest, which can also be 
measured in experiments with ultracold fcrmions 32|, is 
the double occupancy (n-j-n^). At half filling, (n^h^) is 
expected to decrease with decreasing temperature. This 
can be seen in Fig.QJb), where we plot DQMC (symbols) 
and NLCE (lines) results for the double occupancy vs T 
for three values of U in the honeycomb and square lat- 
tices. (Note the excellent agreement between the results 
obtained utilizing the two approaches.) At high tempera- 
tures, {nfhi) is essentially the same for both geometries. 
However, as the double occupancy decreases when re- 
ducing T, one can see that the results in the honeycomb 
lattice depart from, and remain at higher values than, 
those in the square lattice. As this occurs, an upturn 
can be seen in the double occupancy with decreasing T. 
Especially for small U /w, this upturn is more pronounced 
in the honeycomb lattice than in the square lattice (note 



\ (a) 1.04 



— T/w=0.099 (HC) 

— T/w=0.061 (HQ 



U/w=3/2 

T/w=0.102 (SQ) 
T/w=0.063 (SQ) 



0.35 0.55 0.75 0.95 1.15 
\l/w 




3 



FIG. 1. (Color online) (a) NLCE results for the compress- 
ibility vs U in the honeycomb (HC) and square (SQ) lattices, 
at half filling, for two values of T/w that are very close in 
both lattices. Inset- Equation of state for U/w = 3/2 and 
the same two values of T/w as in the main panel. These 
results were obtained after three cycles of improvement of 
Wynn's resummation algorithm [2(|. The zero chemical po- 
tential, which corresponds to half filling for the particle-hole 
symmetric Hamiltonian, is shifted by U /2 for the nonsymmet- 
ric representation of the Hamiltonian in Eq. (TTJ). (b) DQMC 
(symbols) and NLCE (lines) results for (n^h{) vs T in both 
lattices at half filling for U/w = 1/2, 1, and 3/2. Hexagons 
(Squares) and solid (dashed) lines correspond to honeycomb 
(square) lattice. Statistical error bars for DQMC data are 
shown only when they are greater than the symbol size. The 
NLCE results were obtained using Euler resummation, and 
we report the last order (thick lines) and the one to last order 
(black thin lines). 



that for U/w = 1/2, it is absent in the latter geometry). 

The existence of a region in temperature in which there 
is an anomalous d{fi^n{) / dT < has been discussed in 
the context of the Hubbard model in the square lattice. 
Early dynamical mean-field theory calculations identified 
a significant anomalous regi on J33l| , which was later found 
to be marginal in DQMC [Ijand NLCE calcula- 
tions in two dimensions (2D). Interest in the existence 
of such a region developed as it signals adiabatic cool- 
ing with increasing U . This follows from the relation 
dS/dU = —d{n^h^)/dT [33|], which implies that at con- 
stant T, the entropy (S) increases (or, that at constant 
S, the temperature decreases) with increasing U. DQMC 
34 1 and NLCE [29[ calculations have also shown that, 
starting with short-range spin correlation for small values 
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FIG. 2. (Color online) Isentropic curves for the temperature 
vs U /w in the honeycomb lattice at half filling at several con- 
stant entropies. The crossover scale, T* , from NLCE is only 
depicted in the regime where the ground state is a Mott in- 
sulator with long-range AF correlations. For ,5=0.20, we also 
show results for the square lattice (dashed line and square 
symbols) . The inset shows the entropy per particle at T = T* 
vs U/w. Lines (Symbols) correspond to NLCE (DQMC) re- 
sults. 



of U, one can generate exponentially large AF correla- 
tions by adiabatically increasing U , despite the fact that 
there is almost no cooling for weak interactions. How- 
ever, the entropy per particle needs to be S < 0.5. 

We plot in Fig. [2l the isentropic curves in the T — U 
plane for the honeycomb lattice. By comparing those re- 
sults with the ones for the square lattice (see Refs. 2i| 34| 
and the results for S = 0.2 in Fig. [2]), it becomes apparent 



that adiabatic cooling is more significant in the honey- 
comb lattice for small values of U. This occurs in the 
absence of a Mott insulating ground state and where the 
available number of states at any given T in the honey- 
comb lattice is smaller than in the square lattice [see the 
compressibilities in Fig. HJa)]. One may wonder if this 
could ease the realization of exponentially large AF cor- 
relations in the honeycomb lattice in comparison to the 
square lattice, where it remains a major experimental 
goal [35)]. The region with exponentially large correla- 
tions can be identified from T*, which is the tempera- 
ture at which the uniform susceptibility is maximal for 
U beyond the critical value for the formation of the Mott 
insulator. T* is also plotted in Fig. [2] and shows that 
an entropy per particle S < 0.6 is needed to generate 
exponentially large correlations in the honeycomb lat- 
tice. This is close to, but above, the entropy required in 
the square lattice. The entropy per particle at T* in the 
square and honeycomb lattices for half filled systems, S* , 
is shown in the inset of Fig. [2] Beyond U/w = 1 , S* can 
be seen to be almost the same in both lattice geometries 
(S* < 0.5). 

Probing long-range AF correlations turns out to be 
very challenging in optical lattice experiments. As a first 



FIG. 3. (Color online) (main panel) Nearest-neighbor spin 
correlations and (inset) next-nearest-neighbor spin correla- 
tions in the honeycomb and square lattices at half filling as a 
function of entropy, for U = uu/2 and 3w/2. Lines and sym- 
bols are the same as in Fig. Hlb). Note that, for S" nn in the 
square lattice, only QMC results are shown. 



step towards this goal, and towards identifying the AF 
Mott insulator in the Hubbard model on the square lat- 
tice, experiments have already measured NN spin corre- 
3(1, 37 1 . They increase as the temperature 



lations 

is lowered and can be significant even before long-range 
order sets in the system. In Fig. [31 we plot S^ z n in the 
square and honeycomb lattices vs S for two different val- 
ues of U/w, one below and one above the critical value 
for the formation of the Mott insulator in the honeycomb 
lattice. That figure shows that, unexpectedly, there is an 
extended region in entropies where \S^\ are greater in 
the honeycomb lattice than in the square lattice, and that 
this happens even when the ground state in the former 
is a semimetal or a spin liquid while it is an AF Mott in- 
sulator in the latter. At very low entropies, we find that 
, ultimately, \S*„\ in the square lattice becomes greater 
than in the honeycomb lattice, but the entropy at which 
this occurs becomes smaller as U increases. 

Our results imply that strong NN spin correlations can 
be more easily observed in experiments in the honeycomb 
lattice than in the square lattice. They also make evident 
that an enhancement of \S^ z n \ should not be taken as a 
signature of the Neel state, which does not exist in the 
honeycomb lattice for U/w < 0.72, where \S%^\ is greater 
than in the square lattice (for entropies per particle that 
are currently achievable experimentally) . This is because 
such an enhancement can be a very local effect. We have 
also calculated next-nearest-ncighbor correlations, S^ nn , 
in both lattices (see the inset of Fig. [3]) and found them 
to be always stronger in the square lattice than in the 
honeycomb lattice. 

Cooling fcrmions in optical lattices to realize the Neel 
state is currently one of the main experimental chal- 
lenges [35|. To that purpose, one can take advantage 
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FIG. 4. (Color online) (a) Density, (b) entropy, (c) NN 
spin correlations, and (d) compressibility for trapped hon- 
eycomb lattice systems with N — 6.7 x 10 3 particles, an 
average entropy per particle S = 0.67, U = 3w/2, and 
V/w = 1.5 x 10" 3 ,1.0 x 1(T 3 , and 3.8 x 1CT 4 (in order of 
decreasing temperature). DQMC results are depicted as sym- 
bols and NLCE results as lines. Note that, due to the finite 
T grids used in DQMC and NLCE calculations, the values of 
T in both approaches are very close but not identical. 



of the fact that the system is inhomogeneous [a term 



y. Vr 2 

/—iter v i 



where V is the strength of the trapping 



potential and rj is the distance of each lattice site to the 
center of the trap, needs to be added to Eq. (JTJ] and that 
this implies that the entropy is unevenly distributed in 
the gas [3^| . Based on that idea, two recent works, one 
on the square lattice [2!| and the other on the cubic lat- 
tice [HI, have shown that starting with a system with 
high density in the center of the trap (n ~ 2) and with 
an average entropy per particle larger than S* , one can 
achieve a Mott insulator in the center of the trap with a 
local entropy smaller than or equal to S* by adiabatically 
decreasing the confining potential. The excess entropy is 
then stored in the compressible domains with n < 1. 

In Fig. 31 we use the local density approximation 
(LDA), combined with DQMC and NLCE calculations of 
the homogeneous system, to show how the cooling mech- 
anism discussed above works in the honeycomb lattice. 
(For temperatures like the ones studied here, DQMC cal- 
culations have shown that LDA is a good approximation 
on the square lattice Figure 2] depicts the evolu- 

tion of the local density (a), local entropy (b), NN spin 
correlations (c), and the local compressibility (d) as one 
reduces the trapping potential adiabatically in a system 
with U/w = 3/2 and in which the average entropy per 
particle is S = 0.67. This entropy per particle is higher 
than S* = 0.47 for U/w = 3/2. One can see that, as V 
is reduced, the density in the center of the trap changes 
from nearly that of a band insulator to that of a Mott 
insulator [Fig. BJa)]. At the same time, the entropy in 
the Mott insulating region becomes of the order of, or 
smaller than, S* , with the excess entropy being moved to 



the metallic wings [Fig.fJJb)]. This results in strong NN 
correlations in the Mott insulating domain [Fig.[4|c)] and 
a vanishing compressibility in the same region [Fig.^d)]. 
Our results for a specific trapping potential and number 
of particles (similar to the ones in current experiments) 
can be extended to other values of the trapping potential 
and number of particles through the use of the charac- 
teristic density 4fl 41|. 

In summary, we have used DQMC and NLCEs to study 
experimentally relevant thermodynamic properties and 
spin correlations of the Hubbard model in the honey- 
comb lattice. We find that, at half filling and weak inter- 
actions, the compressibility in this lattice may be smaller 
than in the square lattice at low T, despite the fact that 
the ground state in the former is a semimetal and in the 
latter an insulator. We also find that the honeycomb 
lattice exhibits a more significant anomalous region with 
d{h^hi)/dT < than the square lattice, which leads to 
a stronger adiabatic cooling in the former lattice geome- 
try. Remarkably, NN spin correlations in the honeycomb 
lattice are stronger than in the square lattice in an ex- 
tended region of entropies for all U. We discussed how 
these findings are reflected in optical lattice experiments. 
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